tab.path <- "results/table/"
tab.name1<- "app_tab_g07_core.tex"
tab.name2<- "app_tab_g07_end.tex"

yvars <- c("par.onl.sh")
set.seed(2025)

## ParlDiv Clust 
Re$klust <- sapply(strsplit(Re$parl_div, " or "), `[`, 1)
Re$klust <- ifelse(is.na(Re$klust), paste0(Re$county, "-Resid"), Re$klust)


fcontrols <- paste0("dpost:",controls)

df_wb <- data.frame(treat = Re$treat_10, 
  dpost = Re$dpost, 
  yvar = Re$par.onl.sh,
  county = Re$county,
  klust = Re$klust,
  year = Re$year,
  pop =Re$pop,
  dist.roads = Re$dist.roads )

df_wb <- cbind(df_wb, Re[,controls])

for (var in c(controls) ){
  print(var)
  name <- paste0("dp", var)
  df_wb[,name] <- Re[,var]*Re[,"dpost"]
}

df_wb$dptreat <- df_wb$treat * df_wb$dpost
vars <- c( "treat","dpost" ,"dptreat")
fcontrols <- paste0("dp",controls)
df_wb$wry <- df_wb$county == "wr yorkshire"
df_wb$glo <- df_wb$county == "gloucestershire"
df_wb$sur <- df_wb$county == "surrey"
df_wb$nfl <- df_wb$county == "norfolk"

counties <- c("wry", "glo", "sur", "nfl") #Have to create FEs manually for boot_lm command


vars <- c( "dptreat", "dpost", "treat")
outvars <- c( "dptreat", "dpost", "treat")
labvars <- c("DPost X March","DPost", "March")
fcontrols <- c(controls, paste0("dp",controls))


m1 <- felm(as.formula(RegFor( y = "yvar" ,
 x = c(vars,  counties) ,
      FE = "0" , IV="0", clust = "0" )),
           data= df_wb)

boot_lm1 <- boottest(
  m1,  clustid = "klust",
  param = "dptreat",  B = 9999
)
boot_lm1_treat <- boottest(
  m1,  clustid = "klust",
  param = "treatTRUE",  B = 9999
)
boot_lm1_post <- boottest(
  m1,  clustid = "klust",
  param = "dpostTRUE",  B = 9999
)
m2 <- felm(as.formula(RegFor( y = "yvar" ,
 x = c(vars, fcontrols, controls, counties) ,
      FE = "0" , IV="0", clust = "0" )),
           data= df_wb)

boot_lm2 <- boottest(
  m2,  clustid = "klust",
  param = "dptreat", B = 9999
)
boot_lm2_treat <- boottest(
  m2,  clustid = "klust",
  param = "treatTRUE",  B = 9999
)
boot_lm2_post <- boottest(
  m2,  clustid = "klust",
  param = "dpostTRUE",  B = 9999
)
m3 <- felm(as.formula(RegFor( y = "yvar" , 
x = c(vars, controls,fcontrols,counties) ,
      FE = "0" , IV="0", clust = "0" )),
           data= subset(df_wb, year %in% c(1912, 1914)))

boot_lm3 <- boottest(
  m3,  clustid = "klust",
  param = "dptreat", B = 9999
) 
boot_lm3_treat <- boottest(
  m3,  clustid = "klust",
  param = "treatTRUE",  B = 9999
)
boot_lm3_post <- boottest(
  m3,  clustid = "klust",
  param = "dpostTRUE",  B = 9999
)
m4 <- felm(as.formula(RegFor( y = "yvar" , 
x = c(vars, fcontrols, controls, counties) ,
      FE = "0" , IV="0", clust = "0" )),
           data= subset(df_wb, year %in% c(1911, 1912, 1914)))

boot_lm4 <- boottest(
  m4,  clustid = "klust",
  param = "dptreat", B = 9999
)  
boot_lm4_treat <- boottest(
  m4,  clustid = "klust",
  param = "treatTRUE",  B = 9999
)
boot_lm4_post <- boottest(
  m4,  clustid = "klust",
  param = "dpostTRUE",  B = 9999
)
m5 <- felm(as.formula(RegFor( y = "yvar" ,
x = c(vars, fcontrols, controls, counties) ,
      FE = "0" , IV="0", clust = "0" )),
           data= subset(df_wb, pop < 15000 & 
           (dist.roads <= 2 | treat == 1) ))

boot_lm5 <- boottest(
  m5,  clustid = "klust",
  param = "dptreat", B = 9999
)   
boot_lm5_treat <- boottest(
  m5,  clustid = "klust",
  param = "treatTRUE",  B = 9999
)
boot_lm5_post <- boottest(
  m5,  clustid = "klust",
  param = "dpostTRUE",  B = 9999
)

## Spit stars function

SpitStars <- function(x){
  s1 <- ifelse((x <=0.1099), "*","")
  s2 <- ifelse((x <=0.0599), "*","")
  s3 <- ifelse((x <=0.01099), "*","")
  out <- paste0(s1,s2,s3)
  return(out)
}
## Main coef
bsdp <- round(c(coefficients(m1)["dptreat"],
    coefficients(m2)["dptreat"],
    coefficients(m3)["dptreat"],
    coefficients(m4)["dptreat"],
    coefficients(m5)["dptreat"]
), 3)

tstat <- round(c(boot_lm1$t_stat, boot_lm2$t_stat, boot_lm3$t_stat, boot_lm4$t_stat, boot_lm5$t_stat), 3)
pval <- round(c(boot_lm1$p_val, boot_lm2$p_val, boot_lm3$p_val, boot_lm4$p_val, boot_lm5$p_val), 3)
st <- SpitStars(pval)
outcoef <- paste0("$", bsdp, "^{",st, "}$")
outtstat <- paste0("[", tstat, "]")
outpval <- paste0("[", pval, "]")

## Treat
bsdp_t <- round(c(coefficients(m1)["treatTRUE"],
    coefficients(m2)["treatTRUE"],
    coefficients(m3)["treatTRUE"],
    coefficients(m4)["treatTRUE"],
    coefficients(m5)["treatTRUE"]
), 3)
tstat_t <- round(c(boot_lm1_treat$t_stat, boot_lm2_treat$t_stat, boot_lm3_treat$t_stat, boot_lm4_treat$t_stat, boot_lm5_treat$t_stat), 3)
pval_t <- round(c(boot_lm1_treat$p_val, boot_lm2_treat$p_val, boot_lm3_treat$p_val, boot_lm4_treat$p_val, boot_lm5_treat$p_val), 2)
st_t <- SpitStars(pval_t)
outcoef_t <- paste0("$", bsdp_t, "^{",st_t, "}$")
outtstat_t <- paste0("[", tstat_t, "]")
outpval_t <- paste0("[", pval_t, "]")

## Dpost
bsdp_p <- round(c(coefficients(m1)["dpostTRUE"],
    coefficients(m2)["dpostTRUE"],
    coefficients(m3)["dpostTRUE"],
    coefficients(m4)["dpostTRUE"],
    coefficients(m5)["dpostTRUE"]
), 3)
tstat_p <- round(c(boot_lm1_post$t_stat, boot_lm2_post$t_stat, boot_lm3_post$t_stat, boot_lm4_post$t_stat, boot_lm5_post$t_stat), 3)
pval_p <- round(c(boot_lm1_post$p_val, boot_lm2_post$p_val, boot_lm3_post$p_val, boot_lm4_post$p_val, boot_lm5_post$p_val), 2)
st_p <- SpitStars(pval_p)
outcoef_p <- paste0("$", bsdp_p, "^{",st_p, "}$")
outtstat_p <- paste0("[", tstat_p, "]")
outpval_p <- paste0("[", pval_p, "]")

##Output core table
rm(KoreTab)
KoreTab <- data.frame(rbind(
  c("DPost X March", outcoef),
  c("", outpval), 
  vector(mode="character", length=6), #Empty line
  c("DPost", outcoef_p),
  c("", outpval_p),  
  vector(mode="character", length=6), #Empty line
  c("March", outcoef_t),
  c("", outpval_t) 
))

write.table(KoreTab,
    file=paste0(tab.path, tab.name1),
    sep=" & ", 
    quote =  FALSE, 
    col.names = FALSE ,
    row.names = FALSE, 
    eol = " \\\\ \n ")

## Output end of table

MSd <- function(variable){
  x1 <- round(mean(variable, na.rm = TRUE), 2)
  x2<-  round(sd(variable, na.rm = TRUE), 2)
  return(c(x1,x2))
}

Stats <- rbind(MSd(df_wb[,"yvar"]),
  MSd(df_wb[,"yvar"]) ,
  MSd(df_wb[Re$year %in% c(1912, 1914),"yvar"]),
  MSd(df_wb[Re$year %in% c(1911, 1912, 1914),"yvar"]),
  MSd(df_wb[Re$year %in% c(1911, 1912, 1914) & Re$pop < 10000 ,"yvar"])
)
x <- c(
  summary(m1)$r.squared,
  summary(m2)$r.squared,
  summary(m3)$r.squared,
  summary(m4)$r.squared,
  summary(m5)$r.squared)
Rs <- round(x,3)

Ns <- c(m1$N, m2$N, m3$N, m4$N, m5$N)

####OUT
EndTab<- rbind(
AddLines( 6, "Controls", c(LMc("No"),  rep(LMc("Yes"),4))),
AddLines( 6, "Incl. 1913", c(rep(LMc("Yes"),2), c(LMc("No"), LMc("No"), LMc("Yes")))),
AddLines( 6, "Incl. 1911", c(rep(LMc("Yes"),2), c(LMc("No"), LMc("Yes"), LMc("Yes")))),
AddLines( 6, "Pop under 15k", c(rep(LMc("No"),4), LMc("Yes"))),
AddLines( 6, "Within 2~km of roads ", c(rep(LMc("No"),4), LMc("Yes"))),
AddLines( 6, "Mean dep. var." , Stats[,1] ),
AddLines( 6, "Sd dep. var." , Stats[,2] ),
c("Observations", Ns), 
c("$R^2$", Rs)
)
write.table(EndTab,
    file=paste0(tab.path, tab.name2),
    sep=" & ", 
    quote =  FALSE, 
    col.names = FALSE ,
    row.names = FALSE, 
    eol = " \\\\ \n ")

rm(list=
    setdiff(ls()[sapply(ls(), function(x) any(class(get(x)) == 'data.frame'))],
    c("Re","Ind")))
